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Abstract 

Analytical solutions to heat or diffusion type equations are numerous, but there are rather 
few explicit solutions for conditions where the thermal conductivity or diffusion tensors 
are anisotropic. Such solutions have some use in making predictions for idealization of 
real systems, but are perhaps most useful for providing benchmark solutions which can be 
used to validate numerical codes. In this short paper, we present the transient solution to 
the diffusion equation in a cube under conditions of orthotropic anisotropy in the effective 
thermal conductivity or diffusion tensor. In particular, we consider the physically-relevant 
case of transport in a cube with no-flux boundaries for several initial conditions including: 
(1) a delta function, (2) a truncated Gaussian function, (3) a step function, and (4) a 
planar function. The potential relevance for each of these initial conditions in the context of 
validating numerical codes is discussed. 

Keywords: , diffusion, thermal transport, analytical solutions, heat equation, diffusion 
equation 


1. Introduction 

Diffusive-like processes in anisotropic media has many applications, ranging from heat 
transport in composite materials [1] to mass and momentum transport in tissues [2], As 
a recent topical example, magnetic resonance imaging technologies have been developed to 
specifically capitalize on the anisotropic structure of neurological tissues; this has lead to 
a new area of technology known as diffusion tensor imaging (DTI) that has allowed the 
imaging of neurological tissues at unprecedented levels of detail [3, 4, 5, 6, 7, 8]. 

Despite the growing interest in diffusive-like processes in anisotropic media, there are very 
few explicit closed-form solutions available, and fewer that are well-suited for the purposes 
of testing numerical codes. The frequently referenced compendium of Carslaw and Jaeger 
[9, §1.17-1.20] discusses the problem of anisotropy at length, but the only explicit solution 
presented for a rectangular parallelepiped with constant temperature initial conditions and 
zero temperature boundary conditions [9, §14.5]. Even then, the solution is left in integral 


* Corresponding author. 

**e-mail address:brian.wood@oregonstate.edu 


Preprint submitted to International Journal of Heat and Mass Transfer 


July 9, 2015 





form, and the coordinate transformation required to achieve a solution is left to the reader 
to compute. 

In this work, we present explicit series solutions for anisotropic heat transport or diffusion 
in a spatially homogeneous cube for three different initial condition cases. The initial condi¬ 
tions have been selected with the particular purpose of being conditions that are interesting 
from the perspective of validating numerical codes. Because heat transport and diffusion in 
a homogeneous anisotropic medium are mathematically identical, the solutions developed in 
this paper apply to either case. To be specific, note that the diffusion and heat equations 
are given by 


dc Al 


pCp 


dt 

dT Al 

dt 


V ■ (D 7 ■ Vc A7 ) 

(1) 

v • (k 7 • VTA 7 ) 

(2) 


For the case of homogeneous (but anisotropic) systems, we clearly have a correspondence 
between the diffusion tensor, D 7 , and the thermal conductivity tensor divided by the product 
of the density and the heat capacity, K 7 /(pc p ). In the remainder of this paper, we will 
present the solutions for the diffusion equation, with the idea that all of the results are just 
as applicable to heat transport using the simple correspondence discussed. 

In this work, we consider the problem of diffusion of a dilute chemical species in a 
homogeneous but anisotropic medium, whose diffusion tensor can be expressed by a diagonal 
tensor with non-equal entries (i.e., an orthotropic tensor) . Our intent is to provide new 
solutions that focus on novel initial value functions; to date, most of the closed-form solutions 
have been boundary driven. Thus, we have developed solutions for initial conditions that are 
(1) not constant, and (2) of interest for validation of numerical codes. The novel features of 
the material presented are primarily the initial conditions covered, and the explicit, closed 
form of the solutions. 


2. Background 

For the case of diffusion in isotropic systems, a large number of solutions covering many 
possible geometries, boundary condition, and initial conditions can be found in the classic 
texts by Carslaw and Jaeger [9] and Crank [10]. Excellent reviews of more modern papers 
on this topic can be found in references [11] and [12]. For a particular set of boundary and 
initial conditions, we follow explicit solutions as those which contain only sums and products 
of well-known solutions (including conventional algebraic, elementary transcendental, and 
conventional transcendental special functions); following Olver [13], such solutions include 
infinite series. Series solutions can be very useful when they converge reasonably rapidly 
because they provide a method to compute high-accuracy approximations to the solution 
with a finite number of simple computations. Solutions that contain one or more unevaluated 
integrals are terms integral solutions. Although these solutions are also useful, it may not be 
possible to evaluate the integrals that are encountered in terms of algebraic, transcendental, 
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or special functions. Such solutions, ultimately, would require a numerical quadrature to 
evaluate. 

Although the literature on the anisotropic case is reasonably well developed, this litera¬ 
ture has focussed primarily (although not exclusively) on (1) systems that are orthrotropic 
(i.e., the diffusion tensor is diagonal, but the diagonal entries are not equal), and (2) systems 
that are driven by boundary conditions rather than initial conditions. The literature for both 
steady and transient solutions in both 2- and 3- dimensions is summarized in the following 
paragraphs. 

In 2-dimensions, explicit closed-form steady-state solutions have been generated for both 
the orthotropic case [14, 15], and for the fully anisotropic case [16]. For the 2-dimensional 
transient case, separation of variables has been successful for developing explicit closed-form 
solutions for orthotropic cases [17, 18, 19, 20, Ex. 15-6]. The paper by Wang and Chou [21] 
is unique because it is the only result that allows for general initial conditions; all others 
have been been developed for the zero or (equivalently) constant initial condition. A single 
paper by Hsieh and Ma [22] reports explicit closed-form results for a fully anisotropic system 
which are specific to a two-layer infinite slab. 

In 3-dimensions, there are explicit closed-form solutions for the steady orthotropic case 
given by Tungikar and Ran [23], and Hahn and Ozisik [20, Ex. 15-3, 15-4], Aside from the 
well-known solution for the infinite domain given by [9], the only other fully transient explicit 
closed-form solution we were able to find in the literature was the solution of [24] for a finite 
orthotropic functionally graded plate with zero initial conditions. Although Padovan [25] 
outlines a variant of the separation of variables approach for a fully anisotropic transient 
solution, no explicit solutions are provided, and the results are restricted to particular ge¬ 
ometries. To date, there appears to be no general approach that yields explicit closed-form 
solutions for fully anisotropic tensors in which the principle axes are not perpendicular to 
boundaries of the domain. 

Integral solutions (either using Green’s functions or Laplace transforms) have also been 
developed by several researchers for both the orthotropic and the fully anisotropic cases 
in 2-dimensions [26, 27, 28, 29] and 3-dimensions [30, 31, 32, 33, 34, 35, 36, 37]. These 
solutions are formally explicit, but they are not in a closed form. Thus, they these have 
the disadvantage of requiring the computation of integrals before an explicit solution can be 
produced. 

3. Symmetry and Positivity of the Diffusion and Thermal Conductivity Tensors 

In 1851 George Stokes (of fluid dynamics fame) published a paper arguing that thermal 
transport was tensorial in nature, and that the associated transport tensor was symmetric. 
Since that time, there has been substantial discussion on the requirements for the components 
of a conductivity or diffusion tensor to be physically meaningful (see extensive discussions 
on the topic in references [38, 39, 40, 41, 42, 38, 43].) The diffusion tensor (or, equivalently, 
the thermal conductivity tensor) is a positive definite tensor if one accepts the axioms of 
reciprocity laid out by Onsager [44, 45, 41]. Unlike isotropic tensors, anisotropic tensors 
have the potential to have zero transport in one or more principle directions, and thus can 
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only be assured to be positive semidefinte (rather than positive definite). Assuming that the 
tensor takes the form 


Dij, 7 — 


a d e 
d b f 
e f c 


(3) 


then positive semidehnteness is assured by checking Sylvester’s criterion [46] that all of 
the principal minors are all nonnegative, i.e., a > 0,b > 0,c > 0,ab — d 2 > 0 ,bc — f 2 >0, 
ac—e 2 > 0, and a(bc—f 2 )—d(dc—ef)+e(df—be) > 0 . At least one of these quantities must be 
nonzero to prevent the problem from being the trivial case where no transport occurs. Note, 
however, that it is not necessary that all of the principal minors be nonzero. For example, 
in a composite medium constructed of impermeable matrix and spanning capillaries aligned 
in the horizontal plane, the effective diffusion coefficient in the vertical direction is zero, and 
this would be reflected in the principal minors. 

Every three-dimensional second-order positive semideknite tensor has the following prop¬ 
erties: (1) there exist three nonnegative eigenvalues for the tensor; (2) if the three Eigenvalues 
are distinct and positive, then there exist three unique eigenvectors defining the principal 
directions of the tensor; (3) these three principle directions are mutually orthogonal, (4) 
for the case where one or two of the Eigenvalues is zero, or there are positive but repeated 
Eigenvalues, then there are an infinite number of such mutually orthogonal coordinate sys¬ 
tems, and (5) the diffusion tensor is diagonal (but not necessarily isotropic) when expressed 
in a mutually orthogonal coordinate system constructed from its Eigenvectors. Thus, the 
anisotropic diffusion equation given by Eqs. (5)-(7) can always be expressed in terms of a 
orthotropic (diagonal) tensor, D 7 , when put in the proper mutually-orthogonal coordinate 
system; the tensor takes the form 


D 


*J>7 


D xx 0 0 

0 0 

U"uy 

n n D x 


(4) 


Here D xx , D yy , and D zz represent the diffusion coefficients in the direction of the x—, y—, 
and z—directions, respectively. The other two constants represent the anisotropy in the 
system, and are defined by d yy = D xx /D yy , and d 2 z = D xx /D zz . 

For such conditions, the governing differential equation for anisotropic diffusion is well 
understood. The problem in an unbounded domain in 3-dimensions can be stated by 



^ = v • ( o , Vc„ 7 ) 

(5) 

Maximum Principle 

Constraint 

c(x, t) —* 0 as | x| —» oo 

(6) 

I.C. 

c(x, 0) = <h(x) 

(7) 
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This problem requires that the initial condition have finite total mass (or energy); in the 
following we will assume for convenience that the initial condition has a total mass of unity. 

The problem diffusive-like transport in an infinite medium was introduced first by Stokes 
[47]; full solutions to the problem were introduced by Carslaw and Jaeger, although not in 
explicit form. However, following their discussions in §1.17, §2.2, and §10.2, one can, after a 
bit of algebra, extract the solutions for the whole space of the form 


cu 7 ( x ,£) 


d. 


yy 


dz 


f*x=oo ry =oo nz =oo 


2 \Z^D xx t 2\/ tvD xx t 2 \/x l) xx i Jx=— 


$( 2 /, y', z') exp 


oo J y=—oo J z =—oo 


x exp 


(:v - y'fdl 


yy 


4 Dxxt 


exp 


(■g ~ Z ') 2d lz 

4 Dx X t 


dx' dy dz' 


( (x — x') 2 \ 
V / 

( 8 ) 


For the initial condition y', z') = S(x')S(y')S(z') this gives an explicit non-integral solu¬ 
tion of the form. 


c Aj (x,t) 


d. 


yy 


dz 


2\/ 7tD xx t 2 7rD xx t 2^/irD xx t 


exp 


x 


4 Dxxt 


exp 


y 


2 dl 


yy 


^D xx t 


exp 


y2 d 2 ^ 


4:D xx t 


(9) 


which is the frequently referenced solution for an instantaneous point source presented by 
Carslaw and Jaeger [9, §10.2, page 257]. 

In bounded domains, the variable transformations discussed above can still put the dif¬ 
fusion tensor in orthotropic form; however, the boundary conditions are not, in general, easy 
to handle (e.g., separable) in this coordinate system. Thus, coordinate transformation does 
not lead to tractable solutions for the diffusion problem except when the principal axes of 
the diffusion tensor happen to align with the boundaries of the domain (i.e., the problem is 
orthotropic because the domain of the parallelepiped is constructed so that the sides aligned 
with the tensor principle axes.) 

Because much of the terminology and early work on this topic came originally from 
research on the properties of crystals, it is worthwhile to think about an example from a 
crystal system. A monoclinic crystal is a rhomboid where two of the three planes defining 
the system are orthogonal, and the third is not. The natural coordinate system describing 
the geometry of a monoclinic crystal (the crystallographic coordinate system) is one where 
two of the directional vectors are perpendicular, and the third makes an angle other than 
90°. In the natural coordinate system, the thermal conductivity or diffusion tensor can be 
expressed by the anisotropic tensor [9, §1.17] 


D 


*3 ,7 


a<J0 
d b 0 
0 0c 


( 10 ) 


Although an orthotropic representation of the thermal conductivity or diffusion tensor 
in such a crystal does have an orthotropic representation in some coordinate system, that 
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coordinate system is not aligned, in general, with the natural coordinate system of the 
crystal. Hence, the natural boundaries of the crystal would be complicated functions of 
more than one variable (and, hence, unlikely to be easily separable), making it difficult to 
apply separation of variables to the problem directly. This is not a merely academic issue; 
experimental measurements have been made for diffusion in crystals that show exactly this 
behavior. A paper by Bendani et al. [48] describes a set of measurements of the diffusion 
tensor for a naphthalene crystal, and the result is expressed both as an fully anisotropic tensor 
in the crystallographic coordinate system (i.e., analogous to Eq. (10)), and as an orthotropic 
tensor (analogous to Eq. 4) when put in the principle axes of the diffusion process. 

4. Orthotropic Diffusion in a Cubic Domain 

In the remainder of this work, we focus specifically on systems where the coordinate 
system is organized to yield an is orthotropic diffusion tensor, and the boundary conditions 
are defined on planes perpendicular to the coordinate axes. We further assume that the 
domain is a cube with side L; general rectangular parallelepipeds domains can always be 
transformed to a cubic domain as described in the Appendix. Although a few solutions 
do exist in the literature, these solutions have been, to date, either 2-dimensional, or for 
particularly simple initial and boundary conditions. In this work, we illustrate how the 
solution can be obtained for arbitrary (assuming only that the function is an admissible one 
so that the partial differential equation to be well defined) with no-ffux conditions at the 
boundaries. The no-ffux conditions correspond closely to conditions that are of interest both 
for comparison with experimental or numerical results. For this case, the basic mass balance 
equation and boundary conditions take the form 


B.C. 1 
I.C. 


~df 


V ■ (D 7 ■ Vc4 7 ) 


—n 7 ■ D 7 ■ V CA-y = 0, at all bounaries of cube 
c(x, 0) = <3>(x) 


(11) 

( 12 ) 

(13) 


where L = (L,L,L). Note that in the following, we will assume that D 7 is a diagonal 
(but anisotropic) tensor. We will attempt to find a solution using separation of variables. 
Assume CA-y(x,y, z,t) = T(t)U(x,y, z), where U(x,y,z) = X(x)Y(y)Z(z) . Substituting this 
proposed solution yields 

dT 

— U = T D 7 : VW (14) 

dt 7 v ' 

where here VV U = ^ . Note that because D 7 is diagonal, we have D 7 = 0 for i ^ j. 

Thus, we have 


D 7 : VW 


Dij 


d 2 U 

dxjdxi 


d 2 U n d 2 U , n d 2 U 

lhA + ”!h? + “ nr 


( 15 ) 


where each of D xx , D yy , and D zz are (potentially) unique. 
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4-1■ Separation of Variables 

For this problem, the process of separation of variables can be used efficiently to de¬ 
velop fully transient, 3-dimensional solutions. Although the approach is reasonably straight¬ 
forward, the separation does require some careful handling of the separation constants in¬ 
volved. In any event, the separation of this problem does not appear to be well described in 
the literature. As a result, the process is described in detail in this section. 

Using the functional relationship U(x,y,z ) = X(x)Y(y)Z(z ) and substituting, we find 


XYZ°4- = YZD . ** 
dt 


T + XZD„ 


d 2 Y 


dx 2 ' yy dy 2 

Dividing both sides of this equation by XYZT yields 


d 2 Z 

T + XYD ZZ —T 


(16) 


ldT _D xx d 2 X D yy d 2 Y D zz d 2 Z 

T dt X dx 2 + Y dy 2 + Z dz 2 [ } 

For later convenience, we normalize this equation by D xx . Defining d 2 y = D xx /D yy and 
d 2 zz = D xx /D zz , we have 


1 dT _ 1 d 2 X 1 d 2 Y 1 d 2 Z 
D XX T ~dt ~ XlhV + d yy Y dy 2 + d 2 zz Z dz 2 ^ ’ 

Noting that for this equation, the right-hand side is a function of only X, Y, and Z, and the 
left-hand side is a function of only T, we conclude, as conventional, that both sides must be 
equal to a constant. The same argument can be made for each of the terms on the right- 
hand side independently; each term is independent of all of the others. Thus, we define the 
separation constant A 2 by 


1 dT _ 1 d 2 X | 1 d 2 Y | 1 d 2 Z 

D XX T ~dt ~ + d 2 yY dy 2 + d 2 zz Z dz 2 


we note also that, because each of the terms on the right hand side is independent of all 
others, that each of these terms is also equal to a constant. For convenience, we define these 
three constants by 

A 2 = a 2 + f3 2 + y 2 (20) 

This defines a set of four ordinary differential equations as follows 


^ + A 2 D xx T = 0 
dt 

(21) 

d 2 X „ . . 

+ OL X — 0 

dx z 

(22) 

iv 2 + V 1 - 0 

(23) 

d 2 Z „ „ 

dz* + d ’° 7 Z ~ ° 

(24) 
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The solutions to each of these is straightforward. For the time variable, the solution is a 
decaying exponential of the form 

T(t) = T 0 exp (-A 2 D xx t) (25) 

For the remaining equations, the characteristic equation yields roots of the form ±m, 
±dyy/3i , and ±d zz ^i. Thus, the solutions are trigonometric functions of the form 


X(x) = Alsin(ax) + B cos(ax) (26) 

Y(y) = G sin (d yy /3y) + H cos (d yy /3y) (27) 

Z(Z) = Rs'm(d zz ^fz) + S cos(d zz ^z) (28) 

The use of the no-flux boundary conditions require that these functions obey the relationship 

X'(x) = aAcos(ax ) — a/LB sin (ax) (29) 

Y\y) = d yy /3G cos(d yy (3y) - d yy (3H sm(d yy /3y) (30) 

Z\z) = d zz ~/Rcos(d zz ^z) - d zz ^S sin(<4*7*0 ( 31 ) 


For each of the three boundaries where the origin is included, these relationships indicate 
that A = G = R = 0. Thus the solution is of the form 


X(x) = B cos(ax) (32) 

Y(y) = Hcos(d yy /3y) (33) 

Z(Z) = S cos(d zz 'yz) (34) 

The remaining boundary conditions require the following 

sin(aT) = 0 (35) 

sin (d yy f3L) = 0 (36) 

sin(d zz ^L) = 0 (37) 


and these are met, respectively, by Eigenvalues of the form a = £ir/L, d yy /3 = rrm/L , and 
dzzl — nit/L for £, m, n = 0,1, 2 ... The solution of the set of differential equations, then, is 
given by the linear combination of all possible Eigenfunctions. 

£=oo n =oo m —oo 

c A ^(x,y,z,t) = EE E exp (—A 2 D xx t) B(H m S n cos(£iTj;) cos(m7r^) cos(n7r |-) (38) 

£=0 m =0 n =0 

Recalling the relationship A 2 = a 2 + /3 2 + y 2 , then we also have A 2 = (£ 2 + m 2 /d yy +n 2 /d 2 zz )j^. 
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(39) 


l=oo n= oo m= oo / 2 

CA^x,y,z,t) = YYY ex p (-^ 2 + C + £)^^ 

£=0 m =0 n =0 ' 

x B(H m S n cos(£7r£) cos(m7r£) cos(wr£) 

The constants B^, H m , and .S',, are found in the conventional manner using the initial 
condition and the orthogonally of the trigonometric functions. Specifically, for t = 0 we have 

£=oo n=oo m= oo 

y, z) = B^H m S n cos^nf) cos(m7r£) cos(n7r-|) (40) 

£=0 m=0 n =0 

multiplying both sides of this expression by cos^V^) cos(m / 7T^) cos(nV^) and integrating 
yields 


r*x=L 


' x=0 


n y=L 


, y =o 


rz=L 


&(x, y, z) cos(t' / 7rr) cos(m , 7Ty) cos(n'nj-)dz dy dx = 


' 2=0 


£=oo n=oo m=oo „ X= L „ y= L „ Z= L 

Y Y / / / B t H m S n COs(£7Tf) COs(f 7Tf ) 

£=0 m=0 n= 0 ^ J z=0 

x cos(m7r^) cos(mV|) cos(n7r|;) cos(n'7ij^)dz dy dx (41) 

The right-hand side of this expression is non-zero only for the condition £ = m = rn ', n = 
n', for which we have 



fj/=L /-2=L 


<h(a;, ?/, z) cos(£7t|-) cos(m7r|-) cos(m:f )dz dy dx = 


* ?/=0 J z=0 


f*x=L ry=L pz=L 


BtHmSn 


I x=0 J y =0 J z =0 


cos 2 (£7t^) cos 2 (m7r|;) cos 2 (mrj^)dz dy dx 

(42) 


The trigonometric integral on the right-hand side is easily verified to be L 3 /8 when £, m, n > 
0. In general, the value of the integral is L 3 /(2 JV °), where iV 0 is the number of non-zero 
indexes (i.e., 1V 0 = 0, 1, 2 or 3). Thus, the value of B ( H m S n is found from 


2 n ° 

BiH m S n — 


f*x=L ry=L pz=L 


$(a;, y, z ) cos(£7t£) cos(m7r|;) cos(rntjr)dz dy dx (43) 


>x=0 


'y =o 


'2=0 


Note here that the Fourier coefficients are not necessarily independently determinable. For 
non-separable initial conditions, only the product of the coefficients can be determined; 
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frequently a single symbol is used for this coefficient (e.g., B (mn = BpH m S n ). The solution 
is given by 


1=00 m= oo n= oo 


c Al (x,y,z,t) = Y.l>2J2 ex v (-[ £2 7 ? d - + w + g-£ D xx\t 


£=0 m= 0 n=0 


Be m n cos(t'7r£) cos(m7Tf) cos(n7r|;) 
where the Fourier coefficients are evaluated from 


(44) 


2^0 rx=L ry=L r-z=L 

Btmn = —rzr / / / $(x, y, z) cos(£n f) cosirmr y) cosinn f) dx dy dz (45) 

Ja;=0 jy=0 J z= 0 

For the non-separable initial conditions, it is necessary that the initial condition function be 
well-behaved enough such that a representation by a triple Fourier series is possible. For the 
practical purposes discussed here, it is sufficient that the initial condition has meaning in at 
least a distributional sense. 

There is a useful simplification for functions that are multiplicative in the form <F(a;, y, z ) = 
< F z (a;) < f ) ?; (|/) < F 2 (^) (i.e., separable). For this case, we have 


2 n ° 


(*x=L 


ry=L 


f t z=L 


B e H m S n =— $ x (x)cos(£iTy)dx %(y) cos(m7rf) dy 

B l x = 0 Jy= 0 


z=0 


& z (z) COs(mTy) dz 
(46) 


and this allows us to determine By H m , and S n independently by 


Bp = 


2 N ot 


f t x=L 


' x=0 


Hm 


2 n o m ry=L 


S n = 


n Jy=v 

2 n o n rz=L 


' 2—0 


$ x (a;) cos(f7Ty) dx 
$ y (y) cos(rmTy) dy 
® z (z) cos(wr|)) dz 


(47) 

(48) 

(49) 


where here, N 0 p,N 0m , and N 0n are equal to zero if the associated index (£,m, or n) is zero, 
and equal to 1 otherwise. For the problems considered here, we have imposed the condition 
that the total mass in the system is unity; this means that we can compute the zeroth 
coefficient immediately for the multiplicative initial condition directly by 


B 0 

Ho 

So 


1 

L 

1 

L 

1 

L 



/‘Z=L 

' z=0 
r-z=L 

I z=0 
pz=L 

' z =0 


: (x)dx 


® v {y)dy 


z{z)dz 


1 

L 

1 

L 

1 

L 


(50) 

(51) 

(52) 
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Because the series involved is a cosine series, it is convenient to make the first term 
explicit, and represent the remaining terms as a sum from 1 to oo. For the multiplicative 
initial condition, this yields a particularly simple expansion 


c Aj (x,y,z,t) = 


— D t 


X 


Bq + ^ B t cos(^7r^) exp 

t =1 

m —oo 

H o + ^ hT m cos(m7rexp (^-^f^D xx t 

m— 1 
n=oo 

S 0 + ^ 5' n cos(n7rf)exp (~%rj?D xx t 


n= 1 


(53) 


^.2. Moments 

In applications, it is often desirable to have explicit formulae for the spatial moments 
of the solution. One reason for this is to provide a global (spatially) measure that can be 
compared with either experimental data or with numerical results to provide some sense of 
how much agreement there is with the analytical solution. The zeroth, first, and second 
moments are usually three that are considered in applications. The zeroeth moment (total 
mass) in each case is normalized to unity. The first moment measures the center of mass 
of the solution as a function of time, and in the steady state it should be identically equal 
to x = (L/2, L/2, L/2) regardless of the initial condition. The second centered moment is 
particularly useful for diffusion problems, because it relates to the spread of the solution 
around its first moment. 


Zeroth moment 


m 0 (t) = 


rx=L ry=L nz=L 


CA-y(x, y, z, t)dzdydx 


' x=0 


>y =o 


I z =0 


(54) 


For convenience, the zeroeth moment (total mass) for each case is normalized to unity, so 
that mo = 1. To obtain a different value for the total mass, m' 0 , one need only multiply the 
unit mass solution by the appropriate constant so that c'(x,y,z,t) = m' 0 c(x,y, z,t), where 
<d(x,y,z,t) is the concentration associated with the mass m' 0 . 


There are three first moments for the concentration held, one for each of the three Carte¬ 
sian axes. The first moment is found by evaluating the following integrals. 

First moment 


m x (t) = 


m y (t) = 


m z (t) = 


f*x=L ry=L pz=L 

/ / / 

I x=0 J y =0 J z =0 

rx=L ry=L t*z=L 

/ / / 

> #=0 J y =0 J z =0 

f‘x=L ry=L pz=L 

/ / / 

f x=0 J y =0 J z =0 


xcA'yix, y, z, t)dzdydx 


ycAy^x, Vi z, t)dzdydx 


zcA'yix, y, z, t)dzdydx 


(55) 

(56) 

(57) 
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For the case of multiplicative initial conditions, the first moments can be computed explicitly 
by 


m x (t) 

m y (t) 

m z (t) 


L 2 

Y 

L 2 

Y 

L 2 

Y 


1 . 2 (-! + (-!)') 

~l + Y‘2= b ' - p - exp 

1=1 


l 2 YY (-1 + (-1)™) 

7 + ~2 H m -2-“ eX P 

L 7 r z z ' m z 

m= 1 

i + A'f s „(=l±77) exp 

L n z n z 

n —1 


— n t 

j^2 


9 9 

m 7rj_ / 

~ J2 T 2 J-^XX L 

a yy ^ 


(J2 £2 J-'XX 1 ' 


(58) 

(59) 

(60) 


For initial conditions that are symmetric about the point x = (L/2, L/2, L/2) for each of 
the three coordinate axes, the first moments are by definition m x = rn y = m z = L/2; thus 
they do not change in time. 


The centered second moments are actually represented by a 3 x 3 tensor, of the form 


M xx 

M xy 

M xz 

M yx 

Myy 

M yz 

M zx 

M zy 

M zz 


(61) 


In practice, however, it is generally the diagonal elements of the centered moment tensor are 
of primary interest. The diagonal elements of the centered second moment tensor are given 
by the following. 

Centered second moment 


M xx (t ) 


Myy(t) 


M zz [t ) 


r-x=L ry=L pz=L 

/ / / 

f #=0 J y= 0 J z= 0 

rx=L ry=L rz=L 

/ / / 

I x=0 J y =0 J z=0 

fx=L ry=L pz=L 

/ / / < 

' x=0 J y =0 J z =0 


/x — -^) z CA'y(x, y, z , t)dzdydx 


(;y — ^fCA'yix, y , z, t)dzdydx 


z — k) 1 CA- l (x,y, z,t)dzdydx 


(62) 

(63) 

(64) 


For reference, the off-diagonal elements, A/,, h of the centered second moment tensor are 
defined analogously, e.g., 


px=L ry=L pz=L 

M xy = / / / (x - ^)(y - ^)c Al (x,y,z,t)dzdydx (65) 

J x=0 J y —0 J z =0 

however, the off-diagonal elements will not be considered further here. 

For the case of multiplicative initial conditions, an explicit representation of the second 
moments can be determined by evaluating the integrals above. This leads to the closed-form 
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expressions 


M xx (t) 

Myy(t ) 

M zz (t) 


12 

12 

12 


1 12 ^ „ (1 + (-1)*) 


IZi V—A ^ 

H—^ ^ Be 

TV 2 


1=1 
m =oo 


7 2 


exp 


t%B xx t 


(1 + HD 


n2 


2 22 _> 

~ ^ ^2 eX P ^y^h^xxt 


7T 


ra=l 

n=oo 


1 12 5 ( 1 + (" 1 ) w ) 7 n 2 ^ t 

+ n 2 exp ^ dL L 2 

n=l 


( 66 ) 

(67) 

( 68 ) 


It is clear from these expressions that at arbitrarily large times, these moments all tend 
toward the value L 2 / 12, which is the correct value for the centered moment of a uniform 
cube. 


5. Explicit Series Solutions for Particular Initial Conditions 

In this section, the series solution for the four initial conditions described previously are 
computed. The first (when it is different from the value 7) an d second moments are also 
computed as series. To provide a more intuitive feel for the solutions, isosurface plots of the 
concentration held are presented. For presentation of the plots, we define the diffusive time 
scales in each of the cardinal directions by 


rj -\* 


r£ r* 


M 

D X x 

(I ) 2 

d 2 D 

^yy-^XX 

(I ) 2 

dlzDxx 


(69) 

(70) 

(71) 


These represent, very approximately, the time it takes for diffusive gradients in the each 
direction to move a substantial amount of mass to the boundary. Intuitively, gradients in 
the each direction should approach zero after just a few multiples of the diffusive time scale 
for that direction. For the plots generated in this work, we will use the set of parameters 
defined in Table 1; these parameters represent physically reasonable values for diffusion of 
dilute aqueous species in fluids or gels. For other physical systems (e.g., heat transfer, 
diffusion in solids), these parameters would obviously have different values. However, there 
is substantial utility in providing examples that correspond physically to some system, as 
we have done here. As seen in Table 1, T* is the smallest time scale. Thus, for plotting we 
have adopted the following normalized time scale 


t* 


t 

'J 1 * 
- L X 


(72) 
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Table 1: Parameters used in the example computations. 


Parameter 

Description 

Value 

Units 

L 

Domain Size 

0.01 

[m] 

m 0 

Total mass of diffusing solute 1.0 [fig] 



r — mo 

L 'OC p 3 

Steady-state solute concentration 

1 x 10 6 

[gg/m 3 ] 

M 0O = g 

Steady-state second moment 

8.33 x 10“ 6 

[m 2 ] 

D x x 

Diffusion coefficient, x— direction 

1 x 10“ 9 

[m 2 /s] 

d ly 

Anisotorpy ratio, y— direction 

2 

H 

d'L 

Anisotorpy ratio, z— direction 

4 

H 

rp* 

1 X 

Diffusion time scale, x— direction 

25000 

[s] 

nn* 

1 y 

Diffusion time scale, y— direction 

50000 

[s] 

1 Z 

Diffusion time scale, z— direction 

100000 

[s] 

a 

Size parameter for initial condition, Case 2 

0.5L 

[m] 


Variance for initial condition, Case 3 

0.1L 

[nr] 

Ky 

Parameter for planar initial condition, Case 4 

20 

H 

K z 

Parameter for planar initial condition, Case 4 

40 

H 


To provide additional confidence in the series solutions have been computed, each using 
only the first 20 terms of the series, and compared with the results of a previously verified 
numerical code. The code and convergence is described in detail in §6. 

5.1. Case 1 

For Case 1, the initial condition is given by the delta function. The delta function provides 
a good case for assessing how well codes conserve the dependent variable for challenging 
initial conditions, and are also useful for verifying particle tracking codes [49]. The initial 
condition for this case is specified by 


I.C. c(x,Q) = 8(x-L/2) = 8(x-L/2)8(y-L/2)8(z-L/2) (73) 

The values of the Fourier series constants are found by the conventional arguments using 
the orthogonality principles of trigonometric series. This evaluation leads to 
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B 0 = H 0 = S 0 = - 


B e = — j 5{x — L/2)cos(£n^)dx = —cos(ff), for £ > 0 

2 f L 2 

H m — — 5(y — L/2)cos{mnf-)dy = — cos( n ^ L ), for m > 0 

L Jo L 


S n = — J 8(z — L/2)cos(mrf)dz = —cos[}ff), for n > 0 
Hence, the solution for the concentration field can be given explicitly by 


(74) 

(75) 

(76) 

(77) 




1=00 


X 


X 


1 + 2 ^ cos(^) cos(fb rf) exp (—£ 2 ^D xx t 

i= i 
m=oo 

1 + 2 ^2 cos( r f L ) cos(m7T|) exp (^-ff^D xx t 

m=1 
n =oo 

1 + 2 ^ cos( z f) cos(rarf) exp (~^-^D xx t 


72—1 


Isosurface plots of the normalized concentration appear in Fig. I. 1 
The moments for this solution are found as follows. 


Zeroth moment 


First moment 


nx=L ry=L rz=L 

m 0 (t) = / / / <++(+ y, z, t)dzdydx = 1 

J :r=0 J y=0 J z =0 


injl) = m y{t) = "C,Ui = 77 


(78) 


(79) 


(80) 


The zeroeth and first moments are constant for ail cases considered, except Case 4 (where 
the first moment is transient). Where the moments are constant, they are not plotted. 


^^Mathematica scripts for computing each of the four concentration field solutions and associated moments 
as series will be included as ancillary materials in the final submission of this manuscript. 
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Figure 1: Concentration isosurfaces for the 5-function initial condition (Case 1); concentration is normalized, 
= . (a) Although formally not plottable, the initial condition is represented as a point distribution with 
finite mass, (b) The concentration field at t* = \T* using the first 20 terms of the series solution. 


Centered second moments 


px=L ry=L nz=L 

M xx (t) = / / / (x ~ ^) 2 c Aj (x,y,z,t)dzdydx 

J x=0 J y =0 J z =0 

£=oo 


12 


1 + E^l ex P (-t%D xx t) cos (f) (1 + (-1)') 


£ 2 7T 2 

t=\ 

px=L ry=L nz=L 

M m (t) = / / / {y ~ ^) 2 c Ay (x,y,z,t)dzdydx 

J x=0 J y =0 J z=0 


12 


i+ E J^ exp (“CSM cos (^ (1+( - ir) 


m 2 7r 2 

771=1 

r x=L py=L rz=L 


M zz {t) — III (z - ^Yc Al (x,y,z,t)dzdydx 


' x =0 J y =0 J z =0 


12 


24 


n= 1 


1 + E ex P( cos (?)(! + (- 1 )”) 


(81) 


(82) 


(83) 


The centered second moments for Case 1 are plotted in Fig. 2. For these plots, we have 
normalized the moments by their steady-state value, i.e., 

L 2 

Moo = Urn [M xx (t)\ = lim [M yy {tj] = lim [M zz (t)\ = — (84) 

t —>00 t —>00 t —>00 1 2 

For plotting purposes, we have defined the following normalized variables 
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Figure 2: Centered second moments for the (5-function initial condition (Case 1) via series solution (solid 
lines), and numerical computations (symbols). 


Ml 


M, 


M: 


M 0 

M, 


yy 


yy 


Ml 


M 0 


M n 


(85) 

( 86 ) 

(87) 


5.2. Case 2 

For case 2, the initial condition is given by a uniform concentration in a cube of side 
a centered in the domain. This initial condition poses the typical challenges that discon¬ 
tinuities create in numerical approximation. The solution is also complex enough that it 
satisfies the guidelines specified by Roache [50] for defining good analytical solutions for code 
verification. To define the initial condition, we first define the following functions in terms 
of unit Heaviside step functions. 


=H[x - (f - 1)] - H[x - (f + |)] (88) 

*,(») =Hlv - (I - 1)1 - Hit - (f + !)] (89) 

4>«P)=ff[z-(f-l)]-/f[z-(f + f)] (90) 
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Where H(x) is the unit Heaviside step function. The function, <F(x), is given by Defining 
the function <f>(x) by 

$(x) = <$> x {x)<& y {y)<$> z (z) (91) 

The initial condition for a 3-dimensional step function centered in the cubic domain is now 
specified by 


I.C. c(x, 0) = <F(x) (92) 

The values of the Fourier series constants are found in the conventional manner. 


b 0 = h 0 = s 0 = i 

2 f L 4 /“£+f 

Bi = — &(x, y, z) cos(£nj)dx = — / cos(Dr |)dx = 

L Jo L if -1 


4co S (f) S in(g) 


7T a£ 

r L 


, for £ > 0 


. L i a 

2 f 2 [ 2+2 

H m = - ®(x,y,z)cos(rmrf)dy = — / cos(rmrf)dy = 
L io Wf-§ 


4cos(^) sin(ff) 


7T am 
2 rL 


, /or m > 0 


_ L \ a 

2 f 2 + 2 

S n = — $(x, y, z) cos(mrf)dz = — / cos(mrjr)dz = 

L Jo L if-f 


4cos(^)sin(^) 


7ran 


, /or n > 0 


The explicit series expression for the concentration held is 


Cil7 (x,y,^t)=-^ 


1 + 


t^ 4Lcos(f)sin(g) 

7ra£ 

/=i 


cos(Drf)exp ( -f^D xx t ] 


x 


x \ + 4Lcos( ? ) Sm (^) ^ ^ g ■ 

' 7ram L \ d yy L 

m= 1 

7ran 

n=l 

.ce plots of the normalized 
moments for this initial cc 


_ n) 

2L ' cos(m7r/) 

7ram 

7ran — cos(n7rf) exp [~^D xx t) 


n= 1 

e normalized concentration appear in Fig. 3 
this initial condition are given by the follow 


(93) 


(94) 


(95) 


(96) 


(97) 
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Figure 3: Concentration isosurfaces for the step function initial condition (Case 2); concentration is normal¬ 
ized, (a) The initial condition, (b) The concentration field at t* = | T.'* using the first 20 terms of the 
series solution. 


Zeroth moment 


m Q (t) = 


<*x=L 


' x=0 


ry=L 


'y =o 


rz=L 


CA'yix, y, z, t)dzdycLx = 1 


lz= o 


(98) 


First moment 


m x {t) = m y {t) = m z (t) 


L 

2 


(99) 


As with Case 1, these moments are constant in time. The centered second moments are, 
however, transient. Explicit series solutions for these are given by 


Centered second moments 


M xx {t) — ^2 


^ 


M zz {t) = — 


12 


£=oo 


i h— y 'y 

7r 2 ^ 


4Lcos(f)sin(|f) (l+ (-!)*) 


12 


e=i 

m =oo 




i 2 


exp 


?%D xx t 


1 + ^ V 

7T 2 ^ 


4 L cos sin (ff) (1 + (-l) m ) 


12 


m= 1 

72=00 


7 ram 


m A 


i + -2 Y" 

7T 2 


4Lcos(f)sin(ff) (1 + (-!)") 


exp(-gg^f 


2 ,2 


n=1 


nan 


rr 


exp (-%r j?D xx t 


In Fig. 4, the normalized second moments are plotted as a function of time. 


( 100 ) 

( 101 ) 

( 102 ) 


5.3. Case 3 

For case 3, the initial condition is given by a truncated Gaussian function. This function 
is a Gaussian function, centered at x = y = z = L/2 which has been truncated at the 
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Figure 4: Centered second moments for the step initial condition (Case 2) via series solution (solid lines), 
and numerical computations (symbols). 


boundaries of the domain, and normalized so that the total integrated mass in the domain is 
unity. Because the function is analytically smooth, numerical approximations to derivatives 
of it are well-behaved; however, it is complex enough to provide a challenging validation 
benchmark [50]. This initial condition function is specified explicitly by 


exp 


I.C. c(x, 0) = 


( x -§) 


2 crj 




2<rl/dl z 


^2n(j 2 x ^2nal/d 2 y y^2nal/d 2 Z z erf 


k ^ 

2\/2 cr x J 


erf 


2\72((T X ] dyy~j 


erf 


2s/2(cr x /d zz ) f 

(103) 
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(104) 


B 0 = 
B t = 

H m 


‘0 

, 9? 


■S a = 


r / L 2 +2iTr£al \ 

V 2V2La x )_ 

erf (vfc) 

' f / d 2 y L 2 +2iiTma 2 \ 
' l 2 y/2Ld yy a x ) 


COS 


erf 


dyyL \ 
2\/2cr x I 


3? 


erf 


d 2 z L 2 +2inmo 2 
2V2 Ld zz O x 


erf 


(7) exp (-< 2 Sf) , for i > 0 
cos(!f)exp(-^gf), /orm>0 
cosp?)exp(-|-gf), /orn>0 


(105) 


(106) 


(107) 


Note that here at the limit as a x —> 0, this expression reduces to that for the delta initial 
condition. Performing the integrations required to evaluate the Fourier coefficients, the 
closed-form solution for the concentration field is 


c Ay {x,y,z,t) 

1 + 


1 

u 


1=00 Sft 

t= 1 


f ( L 2 +2inla 2 \ 
en| 2V2La x ) 


(108) 


COS 


X 


X 


m= 00 S)fJ 

1+ £2- 

m= 1 


n=oo Sft 

1 + £ 2 - 

n= 1 


erf (vfc;) 

er f f d 2 yy L 2 + 2 ™ m °x \ 


(f) exp ) cos(&rf) exp (-£ : 


,2 hId £ 

JJ2 ^XX u j 


IV 1 a? j 

rtr L y 008 0?) ex P ('€£?) “7m>r!)exp (-^SA»i 

V 2,\/2o x ) 


pr f ( d yy L 
eTI \2y/2a^ 

r / d 2 zz L 2 +2innal 

V 2 V 2 Ld zz a x / 


r ( d zz L 

CTr l 2V2a x 


COS (™) exp cos('U7rf)exp (-%F^D xx t) 


(109) 


Isosurface plots of the normalized concentration appear in Fig. 5. 
The moments for this initial condition are given by the following. 


Zeroth moment 


First moment 


rx=L ry—L rz=L 

m 0 (t) = / / / CA-yiXji/, z,t)dzdydx = 1 

J x=0 J y =0 J z =0 


m x (t) = m y (t) = m z (t) = — 


( 110 ) 


(111) 
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Figure 5: Concentration isosurfaces for the truncated Gaussian initial condition (Case 3); concentration is 
normalized, (a) The initial condition, (b) The concentration field at t* = \T* using the first 20 terms 
of the series solution. 


As before, the first two moments are constant in time. The centered second moments are 
transient, and specified by 
Centered second moments 


L 2 

M X x{t) = Y2 


1 H— o 


24 e= °° 3? 


erf 


L 2 +2in(.cr 2 
2i/2Lcr x 


7 r 


2 

1 =1 


erf 


2\[2cr x 


COS 


7) 


X exp ( fl ' + ' , —— exp (~e%D XI t 


Myy(t) — 


24 


1=00 Sft 


1 + ^ 


£ 2 

erf ( 


dy y L 2 +2iirmcr2 
y 2C^2‘Ldyy(Jx 


m =1 


02; f I dyyL 


2y/2( 

(i + (-i)”) 


«*(t) 




m z 


ex P ( - r ^D xx t 


Mzzit) = ^ 


1 H- o 


24 n=0 ° 3ft 


erf 


d 2 z L 2 -\-2i7rncr 2 
2y/2 Ld zz 


7 r 


2 

n=l 


erf 


dzzL 


Xex P “CTT 


2y/2, 

2 7T 2 <A\ (1 + (- 1 )”) 


COS 


/ 7r n \ 

V 2 / 


77 / 


exp ( - (7 )) D,.,t 


( 112 ) 


(113) 


(114) 


The normalized centered second moments are plotted for reference in Fig. 6. 
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Figure 6: Centered second moments for the truncated Gaussian initial condition (Case 3) via series solution 
(solid lines), and numerical computations (symbols). 


5-4- Case 4 

The initial condition for Case 4 was selected specifically because it is not symmetric about 
the center of the domain. Thus, the center of mass for the solution moves as a function of 
time. Case 4 would be useful, for example, to provide a solution that would detect canceling 
errors in a code; symmetric initial conditions may miss such errors because of the symmetry 
imposed. This case is also the only example provided in this paper for a non-multiplicative 
form for the initial condition. Other non-multiplicative forms can be handled similarly. To 
begin, we specify an initial condition that is a plane, given by 

2 

I.C. c(x,0) = ———— (x + K y y + K z z) (115) 

Ll ~T rvy I r^z) 

Where k is a parameter that increases the slope of the plane along the z— axis. As with 
the previous initial conditions, for convenience the total mass for this function is normalized 
to unity. Unlike previous cases, the evaluation of the Fourier coefficients are not simplified 
by separability. Coefficients are computed by evaluating the integral given by Eq. (45). 
Explicitly, this is 


-Cimn 


2 n ° 


rx=L ry=L nz=L 


L ‘(2 + a) L3 ; I=0 


'y =o 


'2 = 0 


(x+K y y+K z z) cos(£7r£) cos(m7r^) cos(n7r^) dx dy dz 

(116) 
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where £,m,n = 0,1,2,.... A simple integration shows that Be mn is zero when any two of 
the indexes are greater than zero. Thus, the only non-zero contributions for the coefficients 
are of the form Beo 0 , B 0m0 , and B 00n . This yields 


Be oo — 


2^0+1 


rx=L 


L 5 { 1 + Ky + K~) J x= o 


x cos(T7rf) dx = 


4(-l+(—1)*) 1 


Z 3, for £>0 


(l + Kj , + K z ) £ 2 TT 2 

7t— 1 _ l _ —rA, for £ = 0 

(1 + Ky+Kz) L 3 ’ J 


(117) 


BomO ~ 


2 ^ 0+1 


ry=L 


4(—1+(—l) m ) 1 


L 5 (l + Ky + K z ) Jy=0 


K y y cos(m7r-|) dy = 


— ) (l+Kjrt-Kz) m 2 7r 2 


1 


(1+K a +K; 2 ) L 3 ’ 


L 3 ’ 

for m = 0 


for 77i > 0 


(118) 


-Boon — 


2-Wo+i 


r*z=L 


L 5 ( 1 + Ky + K~) J Z= Q 


K z zcos(mnf ) dz = 


4n(—1+(—l) n ) 1 


(1 + Ky+Kz) n 2 7T 2 


1 


(l + Ky+Kz) L 3 ’ 


L 3 ’ 

/or n 


for n > 0 
= 0 


(119) 

Recall, the general solution for the non-multiplicative case is given by Eq. (44); upon sub¬ 
stitution, this yields 


C A j(x,y,Z,t) =-jTg 


1=00 


1 + 


(1 + Ky + Kj) ^ =1 


exp (-£ 2 ^D xx t 


4(-l+ (-!)*) 


£ 2 7T 2 


COs(fbrf-) 


K„ 


(1 + Ky + K z ) 
v y ' m=l 


771=OO 

ex p 1 

771= 1 

{ dl y L 2U ™ X j 

4(-i + (-in 

m 2 7T 2 

77=00 

E ex p( 

71=1 

_n 2 _zZjj A 
<R Z L 2 J 

4(-l + (-l)") 

o o ^ 


cos(mny) 


cos {tvkj) 


( 120 ) 


Note that here, we combined the constant term outside the sum to be consistent with the 
final form of the previous solutions. As a comment, also note that this problem is the 
superposition of the three analogous problems in 1-dimension. This is to be expected because 
of the linearity of the diffusion equation, the boundary conditions, and the initial condition. 
If the initial condition had not been linear, superposition would not have been possible, 
and the result would have been much more complex. Isosurface plots of the normalized 
concentration appear in Fig. 7. 

The moments for this initial condition are given by the following. 


Zeroth moment 


m 0 (t) = 


f x=L ry=L nz=L 


CA'yfx, y, z, t)dzdydx = 1 


' x=0 


'y =o 


'2 = 0 


( 121 ) 
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Figure 7: Concentration isosurfaces for the truncated Gaussian initial condition (Case 3); concentration is 
normalized, (a) The initial condition, (b) The concentration field at t* = |T* using the first 20 terms 
of the series solution. 

First moments Unlike the previous cases, the first moments for this initial condition are not 
constant. Integrating Eq. (120) directly (via Eqs. (55)-(57)), the first moments are given by 


, , L 
"MU =y 


myit) =- 


, , L 
"MU =y 


1 + 


1 + 


^ £=oo 

jz -7 / ex P 

(1 + Ky + K z ) 


7T 2 

- L 2 


— D t 

t 9 J - y XX L/ 


Kn 


(1 + Ky + K z ) 
v y 7 m= 1 




8(-l + (-l)Q 2 

£ 4 ir 4 

8(—1 + (— l ) m ) 2 

m 4 7T 4 


k. 


i + .. v . e “p 8(_1 + <_1)n)2 


(1 + Ky + K Z ) 


n= 1 


n 4 7T 4 


( 122 ) 

(123) 

(124) 


The first moments in this case are different from the first three cases because the center of 
mass of the initial condition is different from the center of mass of the steady-state condition; 
thus, the first moments evolves in time. A plot of the first moments as a function of time is 
provided in Fig. 8. Note for this plot, we have defined the dimensionless variables 


m. 


— 


* m y 

m y=T 

* m z 

m„ =—- 


(125) 

(126) 

(127) 


For this case, the centered second moments are complicated by the transient first moments 
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Figure 8: First moments for the plane initial condition (Case 4) via series solution (solid lines), and numerical 
computations (symbols). 


f‘X=L ry—L pz=L 

M xx {t) — / / / {x - m x (t)) 2 c A ^(x,y,z : t)dzdydx 

J £=0 J y= 0 J z= 0 

px=L ry=L pz=L 

Myy(t) = / / / (V - m y (t)) 2 c Al (x,y,z,t)dzdydx 

J x=0 J y= 0 J z= 0 

px=L ry=L pz=L 

M zz (t) — / / / ( z - rn z (t)) 2 c Al (x,y,z,t)dzdydx 

J x=0 J y= 0 J z= 0 


(128) 

(129) 

(130) 


Substituting Eqs. (120) and (122)-(122) into Eqs. (128)-(128), the closed-form expressions 
for the centered second moment are 

Centered second moments 


M xx (t ) = — - Lm x (t ) + m 2 x {t) 
i 1=00 
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(131) 
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Figure 9: Centered second moments for the plane initial condition (Case 4) via series solution (solid lines), 
and numerical computations (symbols). 
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As a check on this result, one can observe that as time becomes arbitrarily large, each of m x , 
m y , and m z tend toward the value L/ 2. Thus, each of the second centered moments tend 
the to the value ^ ^ ^ = -A-, which the correct value for the steady-state centered 

second moment. A plot of the centered second moments appear in Fig. 9. 
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6. Computations 


To demonstrate the utility of the analytical solutions, we conducted direct numerical sim¬ 
ulations of diffusion in a cube for a prescribed anisotropic effective diffusivity tensor under 
the four different initial conditions using the finite element method (FEM) as coded in the 
commercial package COMSOL Multiphysics 4.4. This code has been previously verified to 
be a correct FEM representation of the diffusion/heat equations (e.g., [49]), in the sense of 
Roache [51]. For the simulations in this paper, the parameters summarized in Table 1 were 
adopted. The code COMSOL solves the discretized equations using the generalized minimal 
residual method with geometric multigrid preconditioning. Successive over-relaxation was 
used in pre- and post-smoothing. To assure that the numerical results were converged, the 
domain was tessellated using a tetrahedral mesh at three levels of refinement. 


We use the Grid Convergence Index (GCI) as proposed by Roache [52] to estimate the 
discretization uncertainties. The GCI is very closely based on Richardson extrapolation and 
has been developed to serve as an uncertainty and convergence analysis tool in computational 
fluid dynamics when analytical solutions are not available. The basic idea is to obtain the 
solution for a given variable 4> on (at least) three grids with spacings Ai, A 2 , and A 3 , 
for which the refinement ratios r 2 1 = A 2 / Ai, and r 32 = A 3 / A 2 , and the variations 
e 32 (x) = </> 3 (x) — </> 2 (x), and e 2 i(x) = </> 2 (x) — </>i(x) are calculated for values of the fine and 
medium grid solutions (</>i and </> 2 ) interpolated on the coarse grid. A local (i.e. pointwise) 
order of accuracy ^(x) can be obtained by solving the following equations iteratively [53]: 


£*(x) 
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+1, if e 32 (x) e 2 i(x) > 0 
-1, if e 32 (x) e 2 i(x) < 0 


The fine-grid convergence index is then computed using: 


(134) 

(135) 


where e \ 1 = |(^i - 4> 2 ) / <P 1 
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Equation (136) can be used to estimate the discretization uncertainty locally or globally. 
& and GCl 2 f ne acquire single values for a global variable (e.g. total mass). When local 
uncertainty estimates are desired, the GCI can be computed locally using a global average 
of & [53] as demonstrated in Fig. 10. Based on these, the numerical uncertainties in the 
fine-grid solutions for the total mass, m 0 , were 0.62%, 0.56%, and 2.06% for the 5-function, 
Gaussian, and step function cases, respectively (data not shown). We define the percentage 
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Figure 10: Sample discretization uncertainty estimates for M xx for the (5-function initial condition; the other 
three initial conditions exhibited similar local errors. Error bars have been magnified lOOx to facilitate 
inspection. 


local discretization uncertainty for the second moment by 




Gdf ine 


x 100 


(137) 


where in this case GCljl ne is the value obtained for computed time point for applicable 
second moment, M xx , M yy , or, M zz . For all for cases the maximum value of /i 2 for all three 
moments was less than 10%. 
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A. Appendix 

In this appendix, we briefly explain the transformation of variables required for or¬ 
thotropic heat or diffusion problems in a rectangular parallelpiped. Coordinate vectors in 
this system are specified by x = (x, y, z), and the parallelpiped has sides of length L in the 
x— direction, and L y , and L z the the other two perpendicular directions. For this case, the 
problem is specified by 
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Now, define the new variables 

(A.4) 
(A.5) 

(A.6) 

From these definitions, it is clear that coordinate vectors ( x,y,z ) span (0,0,0) to (L,L,L) 
as the coordinate vectors ( x,y,z ) span (0,0,0) to (L,Ly,L z ). 

Applying the chain rule twice, it is easy to determine the relationships 
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Finally, denoting 
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we see that the the final form of the governing equation, boundary conditions, and initial 
conditions takes the form 
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(A. 14) 
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where &(x,y,z) = represents the initial condition in the new coordinate 

system. This problem is identical to that in the main body of the paper, establishing that 
for all problems on rectangular parallelpiped, there is an equivalent problems for a cubic 
domain. 
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